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ABSTRACT 

We explain the axisymmetric gaps seen in recent long-baseline observations of the HL Tau 
protoplanetary disc with the Atacama Large Millimetre/Submillimetre Array (ALMA) as be¬ 
ing due to the different response of gas and dust to embedded planets in protoplanetary discs. 
We perform global, three dimensional dusty smoothed particle hydrodynamics calculations 
of multiple planets embedded in dust/gas discs which successfully reproduce most of the 
structures seen in the ALMA image. We find a best match to the observations using three 
embedded planets with masses of 0.2, 0.27 and 0.55 Mj in the three main gaps observed by 
ALMA, though there remain uncertainties in the exact planet masses from the disc model. 

Key words: protoplanetary discs — planet-disc interactions — dust, extinction — 
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1 INTRODUCTION 

Recent long-baseline observations of the HL Tau protoplanetary 
disc at mm to cm wavelengths using the Atacama Large Millime¬ 
tre/Submillimetre Array (ALMA) revealed a striking series of con¬ 
centric, axisymmetric gaps, thought to indicate the presence of 
massive protoplanetary bodies caught in the act of forming plan¬ 
ets (ALMA Partnership 2015). 

While gap carving by planets embedded in discs has been pre¬ 
dicted theoretically since at least Lin & Papaloizou (1986), the ax- 
isymmetry of the rings in HL Tau is puzzling. First, prominent spi¬ 
ral waves induced by the planet are seen in all theoretical simula¬ 
tions of planets carving gaps in gas discs whenever the planet is 
massive enough to open a gap (e.g. de Val-Borro et al. 2006). Sec¬ 
ond, prior observational estimates of the relative disc mass in HL 
Tau suggested a disc of 0.03-0.14 M© around a 0.5-1.3 M© star 
(e.g. Robitaille et al. 2007; Kwon et al. 2011, 2015; the latter au¬ 
thors estimated « 0.1 M© based on fitting disc models to the pre¬ 
vious best observations at mm wavelengths). Massive discs should 
almost certainly show gravitationally instability and hence promi¬ 
nent and large scale spiral arms (e.g. Lodato & Rice 2004), as seen 
in the theoretical predictions for HL Tau made by Greaves et al. 
(2008). Third, HL Tau is in an early phase of star formation (Class 
I evolving to Class II), emitting a 100-200 km/s highly collimated 
jet (Mundt et al. 1990). The ongoing accretion from the envelope 
means that the disc is likely relatively thick, massive and hot, yet 
the clarity of the gaps are more reminiscent of those found in sim- 
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ulations of thin, cold discs (e.g. our most successful attempt to re¬ 
produce the gaps in gas-only discs required H/R< 0.02). 

This begs the question of whether the gaps seen in HL Tau 
are genuine signatures of planet formation, or due to some other 
phenomena such as clumping (Lyra & Kuchner 2013), magnetic 
or Rossby-wave instabilities (Pinilla et al. 2012) or condensation 
fronts (Zhang et al. 2015). Then again, the observation that the first 
four dark rings are located close to what would be resonances be¬ 
tween planets orbiting in those locations, together with slight ec¬ 
centricity observed in the rings are evidence in favour of carving 
by planets (Wolf et al. 2002; ALMA Partnership 2015). 

Key is that the emission seen at the wavelengths observed by 
ALMA is from dust rather than gas. Dust grains settle and mi¬ 
grate in protoplanetary discs at a rate determined by the Stokes 
number, St, the ratio of the stopping time to the orbital timescale 
(Weidenschilling 1977; Takeuchi & Lin 2002) with the stopping 
time depending on the grain size and gas density (Eq. 1). If mm- 
grains have settled to the mid-plane in HL Tau then the gap-carving 
seen by ALMA could be that induced in the thin dust disc rather 
than the relatively thicker gas disc, since it is easier to carve gaps 
in the dust (Paardekooper & Mellema 2004). The absence of spi¬ 
ral structure may also be explained by the dust. Simulations of 
dust grains near gaps carved by planets (Paardekooper & Mellema 
2004, 2006; Maddison et al. 2007; Fouchetetal. 2007, 2010; 
Ayliffe et al. 2012) have demonstrated that dust grains close to 
St = 1 form axisymmetric rings, with smaller dust grains more 
closely following the gas and larger grains trapped in resonances 
(Ayliffe et al. 2012). 

Whether embedded planets in dusty gas discs can explain the 
axisymmetric gaps in HL Tau is the focus of this Letter. 
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2 METHODS 

2.1 Dust/gas simulations 

We perform 3D global simulations of dust/gas discs with em¬ 
bedded protoplanets using the PHANTOM smoothed particle hy¬ 
drodynamics (SPH) code written by DP (Price & Federrath 2010; 
Lodato & Price 2010; Price 2012; Nixon, King & Price 2013). 

We assume a single grain size per calculation. For grains 
of mm-size and larger we employ the two fluid algorithm de¬ 
scribed in Laibe & Price (2012a, b). For calculations with smaller 
grains we employ our new single fluid algorithm based on the ter¬ 
minal velocity approximation (Laibe & Price 2014; Price & Laibe 
2015). This algorithm has been extensively benchmarked on sim¬ 
ple test problems including waves and shocks in dust/gas mixtures 
(Laibe & Price 2011, 2012a; Price & Laibe 2015), but this is our 
first application to planet formation, allowing us to evolve both 
small and large grains efficiently. 
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Figure 1. Constraint on the disc scale height profile obtained from the gap 
widths observed in the ALMA image (upper limits; blue points) combined 
with the lower limit obtained from the dust temperature as a function of 
radius (red points) along with our best fit T oc R~ 0 - 7 (dashed line). 



2.1.1 Drag prescription 


We assume a physical drag prescription where the drag regime is 
chosen automatically depending on the ratio of the grain size to the 
gas mean free path (Kwok 1975; Paardekooper & Mellema 2006; 
Laibe & Price 2012b). All of the grain sizes used in this paper fall 
into the Epstein regime, meaning that the stopping time is given by 

+ _ PgrainS grain /7T7 

kS - n \ o ’ V D 

pCsf V 8 

where fi,,, a u, is the intrinsic grain density (we assume 1 g/cm 3 ), 
Sgrain is the grain size, p — p g + pd is the total density of the 
mixture, c s is the sound speed, 7 = 1 is the adiabatic index and / 
is a correction factor for supersonic drag given by (Kwok 1975) 
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where Av = Vd — v g is the drift velocity between the dust and gas. 
Hence, the Stokes number St = t s Q in the mid-plane is 
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A caveat to the above drag prescription is that it assumes spherical 
grains. Fractal grain shapes or similar would change the effective 
Stokes number associated with a particular grain size. 


2.1.2 Sink particles 

We place three planets, initially located at 13.2, 32.3 and 68.8 au, in 
the three main gaps observed in HL Tau (ALMA Partnership 2015). 
We represent the protoplanetary bodies and the central star using 
sink particles (Bate, Bonnell & Price 1995) with accretion radii of 
0.25, 0.25 and 0.75 au, respectively. Dust and gas can be accreted 
onto the sinks provided that the material is bound to the sink and 
the velocity divergence is negative. The planets are free to migrate 
due to planet-disc interactions and the viscous disc evolution. We 
found migration to be negligible on the timescales simulated. 


2.2 Initial conditions 

2.2.1 Gas 

We set up the discs in PHANTOM using the standard procedure 
outlined in Lodato & Price (2010). The system comprises a cen¬ 


tral star with mass 1.3M© surrounded by a gaseous disc of 10 6 
SPH particles which extends from Ri n = 1 to R ou t = 120 au. 
We assume a power-law surface density profile given by P>(R) = 
Yii n (R/Ri n )- p , where Ei n is the surface density at the inner edge. 
We assume a radially isothermal equation of state P = cl(R)p, 
where c s = c Sj i n (R/Rin )~ q , such that the gas temperature power 
law is T(R) oc ( R/R in )~ 2q and where c Sj i n is the sound speed 
at Ri n . The aspect ratio in the gas is thus prescribed by the p 
and q indices according to H/R — (H/R)[ n (R/Rin)^~ q - We 
set the SPH «av viscosity parameter to 0.1 giving an effective 
Shakura & Sunyaev (1973) viscosity ass ~ 0.005. 


2.2.2 Constraining the disc model from HL Tau 

While Kwon et al. (2011) fit various disc profiles to their previous 
low-resolution observations of HL Tau to constrain the surface den¬ 
sity and temperature profiles, we use the gaps observed in the disc 
(ALMA Partnership 2015) to directly constrain the aspect ratio. 

In particular, while the gap is deeper in the dust (Fouchet et al. 
2007), the width is the same. Since a gap cannot be opened with a 
width less than the disc scale height the three prominent gaps ob¬ 
served in HL Tau provide an upper limit on if as a function of 
radius. Combining this with the dust temperature profile observed 
by ALMA, which, assuming that the dust and gas are thermally 
coupled, provides a lower limit on the gas temperature, we can di¬ 
rectly constrain the sound speed profile. Figure 1 shows the con¬ 
straint provided by the gap widths together with the lower limit 
and our best fitting power law q = 0.35 (T oc R ~° -7 ) with 
(H/R )in = 0.04. This gives (H/R ) ou t ~ 0.08. 

We choose a shallow surface density profile p = 0.1 for com¬ 
putational efficiency since it avoids expensive computation of many 
particles on the shortest orbits. This reproduces the inner disc in our 
simulations, but gives an outer disc that is over-luminous compared 
to the inner disc in the simulated observations (Fig. 4). While more 
work is needed on the disc model, our main aim in the present paper 
is to reproduce the axisymmetry and size of the gaps. 


2.2.3 Gas disc mass and dust-to-gas ratio 

Hayashi et al. (1993) found a H 2 gas mass of 0.03M© contained 
within 1400 au in HL Tau based on the observed 13 CO emission. 
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Figure 2. Rendered image of gas surface density for a disc containing three 
embedded protoplanets of mass 0.2, 0.27 and 0.55 Mj initially located at 
the same distance as the gaps detected in HLTau, at 13.2, 32.3 and 68.8 au. 


We assume this as our starting point, noting the the mass contained 
with a particular radius depends on the density profile according to 


We assume the dust was initially co-located with the gas but has 
migrated to a radius « 10 times smaller than the initial radius to 
give the 120 au dust disc visible with ALMA. Assuming an initial 
dust-to-gas ratio of 0.01 in the undifferentiated material at 1400 
au, this gives E gas ~ Edust within the smaller radius given our 
assumed p index of 0.1. Hence we assume a dust-to-gas ratio of 
unity within 120 au and set Em such that Mdi SC = 0.03M© at 1400 
au. The total gas mass within 120 au is thus 0.0002 M 0 , but this is 
low mainly due to our assumed shallow surface density profile. 

A steeper surface density profile, e.g. p = 1 would give 
a more realistic gas mass of O.OO3M 0 . While even this may 
seem in conflict with earlier fits to much higher disc masses 
(e.g. Robitaille et al. 2007; Kwon et al. 2011, 2015), these assume 
that the dust is co-located with the gas, whereas we expect set¬ 
tling and radial inward migration to produce a more compact dust 
disc, consistent with observations in other systems (Andrews et al. 
2012; de Gregorio-Monsalvo et al. 2013). Moreover, Tamayo et al. 
(2015) found that such a massive gas disc within 120 au would 
smear out the eccentric features observed in HL Tau since the pre¬ 
cession periods of pericentres would be faster than the timescale 
on which planets carve gaps. Our lower disc mass also ensures that 
the mm grains are close to St = 1 in order to guarantee an ef¬ 
ficient migration towards the pressure maxima in the disc induced 
by the embedded protoplanets. Again, a steeper surface density pro¬ 
file would produce less coupling of the mm grains in the outer disc, 
giving dynamics similar to those we find with cm grains. Indeed, as 
we shall see, the spiral waves induced in larger grains by the outer 
planet do appear to be visible at mm wavelengths in HL Tau. 

2.2.4 Dust 

The dust disc is setup assuming an initially constant dust-to-gas 
ratio (equal to unity, see above). Since we simulate one grain size 
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at a time we assume a power-law grain size distribution given by 

n(s) oc s~ m for Smin < S < Smax 5 (5) 

where we assume s m i n — 1/im, s max — 10cm and m — 3.5 and 
perform six different simulations using s = lpm, 10/rni, 100 pm, 
1mm, 1cm and 10cm, respectively (Fig. 3). The grain size distribu¬ 
tion is appropriate to interstellar dust and particles growing by ag¬ 
glomeration in protoplanetary discs (Draine 2006). While approx¬ 
imate, this is more realistic than assuming a constant dust-to-gas 
ratio for each grain size, capturing the main idea that small grains 
are more abundant than large grains. 

In our two fluid simulations we set up the dust disc as a sep¬ 
arate set of dust particles. We use one quarter of the number of 
dust particles compared to gas to prevent dust from becoming ar¬ 
tificially trapped below the gas resolution (Laibe & Price 2012b; 
Ayliffe et al. 2012). For the calculations employing the single fluid 
algorithm (Price & Laibe 2015) there is only one set of particles de¬ 
scribing the mixture, set up as described above for the gas but with 
a dust fraction e set for each particle to provide an initially constant 
dust-to-gas ratio (i.e. e = pd/p g /[l + pd/pg])- With both meth¬ 
ods the dust is then free to settle and migrate in the disc, and feels 
gravity from both the central star and the embedded protoplanets. 

2.3 Simulated ALMA observations 

We perform simulated observations of our model using the 
RADMC-3D Monte Carlo radiative transfer code (Dullemond 
2012) together with the Common Astronomy Software Applica¬ 
tion (CASA) ALMA simulator (version 4.1), focussing on ALMA 
band 6 (continuum emission at 233 GHz). The source is located in 
the position of HL Tau, adopting the disc inclination and position 
angle given by ALMA Partnership (2015). The procedure is as in 
Dipierro et al. (2015) except that the dust distribution is used di¬ 
rectly from our simulation data rather than being prescribed. The 
dust model consists of spherical silicate grains with optical con¬ 
stant for magnesium-iron grains taken from the Jena database 1 . 
Full-resolution images produced by RADMC-3D simulations are 
used as input sky models to simulate realistic ALMA observations 
accounting for thermal noise from the receivers and the atmosphere 
and assuming a perfect calibration of the visibility measurements. 
A transit duration of 8 hours is used to reach an optimal signal- 
to-noise ratio. The observation parameters are chosen to match the 
spatial resolution of the observations (see caption of Fig. 4). 


3 RESULTS 

Fig. 2 shows the gas surface density, while Fig. 3 shows rendered 
images of the dust surface density from the disc model with three 
embedded protoplanets of masses of 0.2,0.27 and 0.55 Mj, respec¬ 
tively, shown after ~10 orbits of the outer planet. None of the plan¬ 
ets produce clean gaps in the gas disc (Fig. 2), but as expected, the 
dust density perturbations induced by protoplanets depend strongly 
on the grain size (e.g. Fouchet et al. 2007; Ayliffe et al. 2012), or 
more specifically on the Stokes number. The density distribution of 
pm -sized particles (top left and centre) are similar to the gas dis¬ 
tribution due to the stronger coupling. The micrometer sized grains 


1 http://www.astro.uni-jena.de/Laboratory/ 
Database/databases.html 
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Figure 3. Rendered images of dust surface density for a disc containing three embedded protoplanets of mass 0.2, 0.27 and 0.55 Mj initially located at the 
same distance as the gaps detected in HL Tau. Each panel shows the simulation with gas plus grains of a particular size (as indicated). 


(top left) capture the spiral density wave launched by the proto¬ 
planets. Millimetre-sized particles (bottom left) are most affected 
by the density waves induced by the protoplanets, exhibiting the 
largest migration towards the gap edges. Interestingly, the dust den¬ 
sity distributions of 10 mm and cm-sized particles (bottom centre 
and right) show axisymmetric waves launched by the planets prop¬ 
agating across the whole disc. The gaps carved by protoplanets in 
cm grains (St ~ 10) show the formation of horseshoe regions. 
As expected, the formation of the horseshoe region in planetary 
gaps depends on the dust-gas coupling, since poorly coupled large 
particles tend to have frequent close encounters with the planet. 
We stress that the important parameter is the Stokes number rather 
than the actual grain size, so equivalent dynamics can be produced 
in different grains by changing the gas density or with different as¬ 
sumptions about the grain composition (see Eq. 3). 

An intriguing possibility is to infer the disc viscosity from the 
morphology of the gaps. Since our results show that the planets are 
not massive enough to open gaps in the gas or, equivalently, the 
viscosity is high enough to effectively smooth the density profile, 
we can only evaluate a lower limit for ass below which planets 
are able to carve gaps in the disc model adopted here. Using the 
criterion in Crida et al. (2006) we infer ass > 0.0002. 

Fig. 4 compares the ALMA simulated observation of our disc 
model at band 6 (right) with the publicly released image of HL 
Tau (left). The pattern of bright and dark rings is readily detectable 
with ALMA using a selection of observation parameters that ensure 
a resolution close to the one reached in the real ALMA image. As 
expected, the emission probes the mid-plane disc surface density in 
large grains (0.1-10 mm), due to their high opacity at these wave¬ 
length (Dullemond et al. 2007; Williams & Cieza 2011). However, 
the intensity in the ALMA simulated observation of this model is 


lower than in the HL Tau image. This difference is likely due to the 
low dust mass, suggesting that the higher disc mass from a steeper 
surface density profile would produce better signal-to-noise ratio. 


4 SUMMARY 

Inspired by the recent ALMA image of HL Tau, we have performed 
3D dust/gas SPH simulations of a disc hosting multiple planets 
aimed at reproducing the peculiar morphology of alternating dark 
and bright rings shown in the observation. We used 3D radiative 
transfer modelling to produce synthetic ALMA observations of our 
disc model in order to test whether dust features induced by em¬ 
bedded protoplanets can explain those detected. We find: 

(i) The absence of spiral structure observed in HL Tau is due to 
the different response of the dust compared to the gas to embedded 
planets. Axisymmetric rings such as those observed are a natural 
product of embedded planets in the disc when observing the dust. 
Spiral structure may still be present in the gas. Hence dust/gas mod¬ 
elling is crucial for interpreting these systems. 

(ii) In agreement with previous authors, we find it is easier to 
carve gaps in dust than in gas. Hence gaps observed in the contin¬ 
uum at mm-wavelengths do not necessarily indicate gaps in gas. 

(iii) We can reproduce all the major features of HL Tau obser¬ 
vations with three embedded planets at 13.2, 32.3 and 68.8 au with 
planet masses of 0.2, 0.27 and 0.55 Mj, respectively, though there 
is uncertainty in the exact masses from the assumed disc model. 
This supports the conclusion that embedded protoplanets are re¬ 
sponsible. Our masses are similar to the 0.2 Mj suggested for all 
three planets in Dong et al. (2014) and consistent with a 10Mj up¬ 
per limit for the third planet given by Testi et al. (2015). 
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Figure 4. Comparison between the ALMA image of HL Tau (left) with simulated observations of our disc model (right) at band 6 (continuum emission at 233 
GHz). Note that the color bars are different. The white colour in the filled ellipse in the lower left corner indicates the size of the half-power contour of the 
synthesized beam: (left) 0.035 arcsec x 0.022 arcsec, P.A. 11°; (right) 0.032 arcsec x 0.027 arcsec, P.A. 12°. 


(iv) The dust density structures for grains with St > 10 (cm-size 
and larger in our simulations; bottom centre and right panels of Fig. 
3) show axisymmetric features that can be identified as waves ex¬ 
cited by the embedded protoplanets (Goldreich & Tremaine 1980). 
These waves can freely propagate across the dusty disc without be¬ 
ing damped by dissipative phenomena (Rafikov 2002). The lower 
drag makes these waves more visible since the drag effect on the 
radial dust motion is negligible. We therefore suggest that the ax¬ 
isymmetric perturbations in the ALMA images in the outer part of 
the second and third planet orbits are generated by the lower cou¬ 
pling of mm-particles. Further simulations with a steeper surface 
density profile should be able to reproduce this. 

There are many caveats to our results, and these are prelimi¬ 
nary investigations. Nevertheless they illustrate that the ALMA ob¬ 
servations of HL Tau can be understood from the interaction be¬ 
tween dust, gas and planets in the disc. 
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